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PREFACE 

Dynamic  partitioning  is  an  array  language  construct 
which  was  conceived  during  the  design  of  the  OL/2  language 
by  the  author.   It  is  interesting  because  it  allows  one  to 
write  array  algorithms  for  areas  which  have  previously  been 
difficult  or  impossible;  for  instance,  the  direct  methods  of 
linear  algebra.   However,  dynamic  partitioning  is  not  re- 
stricted to  array  languages  with  an  algebraic  structure  like 
OL/2.   The  basic  ideas  and  the  information  structure  are 
applicable  to  other  array  languages  with  different  data  types 
and  different  array  operations.   The  purpose  of  this  paper 
is  to  provide  a  formal  description  of  dynamic  partitioning 
and  to  document  the  basic  algorithms  in  a  language  indepen- 
dent form.   This  should  be  valuable  to  anyone  who  is  con- 
cerned with  the  construction  and  design  array  languages. 

The  first  part  of  this  paper  will  appear  in  CACM  [7], 
while  the  latter  part,  which  is  concerned  with  algorithms 
for  implementing  partitioning,  will  be  described  in  this 
paper  in  a  Knuth-like  notation  [6] . 

Dynamic  partitioning  was  first  conceived  in  the  latter 
part  of  1969;  since  then  several  persons  have  been  instrumen- 
tal in  its  implementation.   The  first  was  John  Gaffney,  who 
made  a  feasibility  test.   The  second  was  H.  C.  Adams,  who 
implemented  the  present  version  for  OL/2  as  a  Master's 
thesis  [1] .   All  of  this  would  have  been  impossible  without 
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the  efforts  of  Robert  Bloemer  and  Dale  Jurich,  who  have  been 
major  contributors  to  the  implementation  of  OL/2.   I  would 
like  to  express  my  appreciation  to  each  of  these  people  for 
their  contributions,  confidence,  and  enduring  spirit! 
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ABSTRACT 

The  classical  process  of  partitioning  an  array  into 
subarrays  is  extended  to  a  more  useful  array  language  opera- 
tion.  Various  modes  of  partitioning  are  defined  for  different 
types  of  arrays,  so  that  subarrays  may  vary  over  the  original 
array  in  a  nearly  arbitrary  manner.   These  definitions  are 
motivated  with  several  realistic  examples  to  illustrate  the 
value  of  partitioning  for  array  languages. 

Of  general  interest  is  the  data  structure  for  partitioning 
This  consists  of  dynamic  tree  structures  which  are  used  to 
derive  and  maintain  the  array  control  information.   These  are 
described  in  sufficient  detail  to  be  of  value  in  the  design 
of  other  array  languages.   The  description  presented  in  this 
paper  is  implemented  in  a  new  array  language,  OL/2,  currently 
under  development  at  the  University  of  Illinois. 
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Introduction 

Partitioning,  in  the  usual  sense,  consists  of  separating 
an  array  into  subarrays  by  inserting  partition  lines  between 
certain  rows  and  columns.   The  purpose  is  to  implicitly  define 
the  subarrays  by  their  relative  position  in  the  original  array 
so  that  they  may  be  renamed  or  indexed  for  notational  purposes. 
In  this  primitive  form  partitioning  is  limited,  but  in  a  more 
general  setting  it  becomes  a  very  powerful  array  operation. 

In  essence,  dynamic  partitioning,  as  defined  here  and  in 
the  OL/2  language  [7],  allows  the  partitioning  lines  to  vary 
over  the  array.   This  in  turn  allows  various  subarrays  to  be 
referenced  in  array  expressions.   If  we  also  include  the  option 
of  selecting  subarrays  which  vary  in  size  and  shape,  then  most 
direct  methods  of  linear  algebra  can  be  implemented  in  a  simple 
and  straightforward  manner.   It  is  also  interesting  to  observe 
that  an  array  language  with  partitioning  can  use  storage  as 
efficiently  as  an  element  language  such  as  PL/1  or  ALGOL.   The 
crucial  point  is  not  the  array  language  but  the  algorithm  itself 
If  the  algorithm  is  amenable  to  array  notation  and  temporary 
arrays  are  not  required  by  the  array  expressions,  then  the  al- 
gorithm may  be  implemented  with  a  properly  designed  array  lan- 
guage which  uses  the  same  amount  of  storage  and  the  same  number 
of  operations  as  an  element  language.   To  support  this  point  we 
present  in  the  appendix  two  well-known  algorithms  using  the 
OL/2  language:   the  Crout  decomposition  algorithm  and  the 
Cholesky  method. 
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The  purpose  of  this  paper  is  to  describe  the  various  modes 
of  partitioning  and  a  strategy  for  implementing  it.   A  detailed 
description  of  the  implementation  for  the  OL/2  language  on  the 
IBM  360/75  appears  in  [1] .   Other  array  languages  usually  have 
some  provision  for  selecting  subarrays.   In  particular,  there  is 
the  "take"  and  "drop"  functions  in  APL  [  5  ] ,  and  the  "ranges"  and 
extraction  operators  ^  and  ^  in  the  article  by  Bayer  and  Witzgall 
[  2 ] •   The  SPEAKEASY  language  [  3  ]  and  the  MPL  language  [  4  ]  are 
examples  of  other  array  languages  which  have  a  simple  form  of  sub- 
array  selection. 

Basic  Data  Types 

The  basic  types  of  arrays  for  which  partitioning  is  defined 
are  given  in  Table  I.   From  this  basic  set  a  larger  class  of 
arrays  may  be  constructed,  such  as  block  tridiagonal  matrices  or 
band  matrices.   Even  though  partitioning  applies  to  these  com- 
posite arrays,  it  is  sufficient  to  consider  a  partitioning  strategy 
based  only  on  the  arrays  in  Table  I. 

The  types  of  arrays  which  have  been  chosen  are  the  common 
forms  that  are  encountered  when  solving  problems  in  linear  algebra. 
Since  the  storage  for  arrays  is  usually  critical,  it  becomes  impor- 
tant not  to  store  the  elements  which  are  theoretically  zero.   The 
arrays  still  operate  on  each  other  according  to  the  rules  of 
linear  algebra,  but  the  storage  and  number  of  operations  is  re- 
duced to  a  minimum.   As  a  matter  of  convenience,  the  elements  of 
all  arrays  are  stored  linearly  in  row  major  order.   This  storage 
scheme  together  with  the  various  types  of  arrays  determine,  in 
part,  the  information  that  must  be  inserted  into  the  data  structure 
which  is  described  in  subsequent  sections. 


TABLE  I  Basic  Array  Types 


ARRAY  TYPE 
T 

VIRTUAL 
BOUNDS 

STORED  ELEMENTS 
S(T) 

I 

J 

STRICTLY  LOWER 
TRIANGULAR 

1:n 

1:n 

1  <  J  <  I  <  N 

LOWER  TRIANGULAR 

l:N 

l:N 

1  <  J  <  I  <  N 

DIAGONAL 

l:N 

l:N 

1  <  I  =  J  <  N 

STRICTLY  UPPER 
TRIANGULAR 

1:N 

1:N 

1  <  I  <  J  <  N 

UPPER  TRIANGULAR 

l:N 

l:N 

1  <  I  <  J  <  N 

TRIDIAGONAL 

1:n 

l:N 

Il-J|  <  U 

1  <  I  /J<  N 

RECTANGULAR 

1:M 

1:n 

1  <  I  <  M; 

1  <  J  <  N 

4 
Modes  of  Partitioning 

The  simplest  form  of  partitioning  is  referred  to  as  static 
partitioning.   It  is  defined  by  fixing  the  partitioning  lines 
permanently  after  specific  rows  or  columns  of  an  array.   Figure  1 
illustrates  this  type  of  partitioning  for  several  of  the  arrays 
in  Table  I. 

More  precisely,  we  define  partitioning  for  all  of  the  tri- 
angular and  diagonal  arrays  of  Table  I  by  having  the  partitioning 
line  reflect  off  of  an  imaginary  diagonal  line  which  starts  in 
the  upper  left  corner  of  the  array  and  proceeds  to  the  lower 
right  corner.   Since  these  arrays  are  square,  by  definition,  the 
effect  is  to  simultaneously  partition  the  array  after  the  same 
row  and  column. 

Rectangular  arrays,  on  the  other  hand,  are  partitioned  in 
the  usual  way.  That  is,  they  are  partitioned  after  any  row  or 
column  of  the  array  with  the  row  and  column  partitioning  lines 
being  independent  of  each  other. 

Another  mode  of  partitioning,  which  is  referred  to  as 
diagonal  partitioning,  is  also  defined.   To  simplify  matters 
assume  that  the  array  to  be  partitioned  is  a  square  matrix  of 
order  n.   Then  one  may  specify  any  part  of  the  matrix  which  is 
compatible  with  one  of  the  triangular  or  diagonal  arrays  in 
Table  I;  for  example,  one  may  specify  the  strictly  upper  tri- 
angular part  of  the  matrix.   In  this  instance,  the  effect  is 
equivalent  to  placing  a  diagonal  partitioning  line  just  above 
the  principal  diagonal  of  the  matrix. 
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FIG.  1.   Elementary  partitioning  of  the  basic  types 
of  arrays:   (a)  lower  triangular;  (b)  tridiagonal  ; 
(c)  diagonal;  (d)  strictly  upper  triangular;  (e)  diag 
onal  partitioning;  (f)  a  combination  of  rectangular 
and  diagonal  partitioning  with  a  selection  of  a 
subset  of  elements  represented  by  the  solid  marks. 
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In  all  cases  of  partitioning,  no  additional  storage  is 
needed.   Partitioning  merely  allows  one  to  select  subarrays  of 
various  sizes  and  shapes.   These  subarrays  are  then  available 
for  referencing  in  array  expressions,  so  that  array  algorithms 
may  be  expressed  in  a  simple  and  direct  manner. 

We  may  extend  the  partitioning  definitions  to  the  dynamic 
realm  by  allowing  the  partitioning  lines  to  be  defined  by  ex- 
pressions.  Thus,  in  dynamic  partitioning  the  subarrays  vary 
over  the  parent  array  in  an  arbitrary  manner  as  the  values  of 
the  expressions  change.   At  the  same  time,  the  subarrays  remain 
uniquely  defined  by  their  relative  position  in  the  parent  array, 
so  that  their  name  remains  unchanged  in  any  array  expression. 

To  illustrate,  consider  the  following  statements  which  are 
extracted  from  one  of  the  examples  in  the  appendix. 

LET  A  BE  A  MATRIX  OF  ORDER  (N) ;  FOR  K  =  1,2,...,N; 
PARTITION  A  AFTER  ROW  K-l  AND  AFTER  COLUMNS  K-l,  K; 
SET  C  =  A<2,1>,  X  =  A<1,2>,  Y  =  A<2,2>,  B  =  A<1,3>, 
AND  M  =  A<2,3>; 

The  matrix  A  is  partitioned  into  six  subarrays  in  Figure  8a 
by  the  PARTITION  statement,  and  five  of  the  subarrays  are  renamed 
in  a  corresponding  SET  statement.   The  subarrays  of  A  are  specif ie 
with  the  aid  of  special  parentheses  <  >  so  that  the  common  paren- 
theses can  be  reserved  for  specifying  the  elements  of  the  matrix. 
Thus,  using  conventional  indexing,  A<1,1>  refers  to  the  subarray 
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in  the  upper  left  corner,  and  A<2,1>  to  tne  subarray  immediately 
below  A<1 ,1>. 

In  the  above  example  the  partitioning  is  dynamic,  since  the 
FOR  statement  alters  the  expressions  which  define  the  position  of 
the  part itioning  lines.   As  the  value  of  K  is  altered,  the  sub- 
arrays  C,  X,  Y,  and  B  move  over  the  matrix  A  in  a  uniform  manner. 
The  meaning  of  the  statements  above  for  the  cases  K  =  1  and  K  =  N 
are  discussed  in  the  last  section,  which  is  concerned  with  the  de- 
finition of  null  operands. 

We  now  extend  the  operation  of  partitioning.   Since  the  par- 
titioning lines  implicitly  define  the  subarrays ,  each  subarray  is 
a  candidate  for  further  partitioning.   In  general,  partitioning 
may  be  applied  to  all  subarrays  which  result  from  previous  par- 
titionings,  whether  they  have  been  explicitly  specified  in  a  SET 
statement  or  not.   The  process  is  referred  to  as  nested  partition- 
ing. 

One  of  the  advantages  of  nested  partitioning  is  that  over- 
lapped subarrays  of  various  sizes  and  shapes  may  be  referenced. 
That  is,  the  programmer  has  equal  access  to  any  of  the  arrays 
which  have  resulted  from  nested  partitioning,  provided  only  that 
the  reference  is  within  the  scope  of  the  PARTITION  statement. 
However,  it  is  interesting  to  note  that  in  practice  most  algorithms 
require  only  a  few  levels  of  partitioning.   For  instance,  the 
algorithms  of  Cholesky  and  Crout,  in  the  appendix,  require  only 
one  and  two  levels  of  nested  partitioning  respectively. 

Finally,  there  is  an  even  more  general  set  of  subarrays  that 
may  be  selected  by  overlaying  several  independent  partitionings . 
This  is  achieved  by  first  defining  two  arrays  or  subarrays  to  be 
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equivalent;  for  instance,  SET  T  =  A.   This  means  that  both  T 
and  A  refer  to  the  same  set  of  elements.   Then,  by  allowing  T 
and  A  to  be  partitioned  independently  of  each  other,  we  obtain 
independent  sets  of  subarrays. 

In  conclusion,  we  have  two  basic  types  of  partitioning: 
static  and  dynamic.   With  either  of  these  types  we  may  have 
nested  partitioning  to  an  arbitrary  level.   At  each  level,  there 
is  a  mode  of  partitioning  corresponding  to  a  classification  of 
array  types.   That  is,  one  mode  for  triangular  and  diagonal 
arrays,  another  mode  for  rectangular  arrays  which  is  the 
classical  rectangular  partitioning,  and  finally  a  diagonal 
partitioning  mode  which  selects  triangular  or  diagonal  subarrays. 
Complementing  these  options  is  the  ability  to  define  one  or  more 
independent  partitions  on  the  same  array  by  overlaying  the  par- 
titions . 

The  Data  Structure  for  Partitioning 

Partitioning  is  implemented  by  constructing  a  general  tree 
structure  T.   The  root  of  the  tree  T  corresponds  to  some  array 
A  which  is  explicitly  declared  by  a  LET  statement,  while  the 
subarrays  of  A  form  the  roots  of  the  subtrees  of  T  at  the  next 
level.   If  any  of  the  subarrays  are  partitioned,  then  nodes  are 
added  to  the  subtree  at  the  next  level.   Therefore,  in  general, 
the  entire  data  structure  consists  of  independent  tree  structures 
one  for  each  explicitly  declared  array. 

In  practice,  it  is  not  necessary  to  build  the  entire  tree 
structure;  in  fact,  it  is  preferable  to  build  only  those  parts 
which  are  going  to  be  referenced  by  the  programmer.   This  can  be 
accomplished  if  an  explicit  specification,  such  as  the  SET 
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statement  in  the  last  section,  is  required  with  each  partition 
statement . 

The  justification  for  this  convention  is  based  on  several 
observations.   First,  in  many  algorithms  the  programmer  needs 
to  reference  only  a  few  of  the  subarrays  which  result  from 
partitioning,  and  secondly,  the  size  of  the  tree  structure  can 
grow  very  rapidly  if  all  of  the  subarrays  are  represented  by 
nodes  in  the  tree.   An  additional  advantage,  which  perhaps  is 
equally  important  from  a  semantic  point  of  view,  is  that  the 
SET  statement  helps  to  document  the  algorithm  in  a  natural  way. 

In  some  cases,  however,  it  is  convenient  to  reference  sub- 
arrays  which  have  not  been  explicitly  specified  in  a  SET  state- 
ment.  In  this  case,  the  tree  structure  must  be  altered  dynamically 
during  program  execution.   The  partitioning  lines  still  determine 
the  number  of  subarrays  and  their  relative  position  in  the  parent 
array,  but  it  is  no  longer  possible  to  determine,  a  priori,  which 
subarray  is  going  to  be  referenced.   For  example,  if  A  is  the 
parent  array,  then  A<I,J>  is  a  reference  to  a  subarray  that  is 
not  known  during  compilation.   However,  it  also  follows  that  the 
size  and  shape  of  the  subarray  A<I,J>  may  now  change  as  the  values 
of  the  subscripts  change.   This  provides  a  rather  powerful  form  of 
referencing  subarrays  when  combined  with  dynamic  partitioning. 

It  is  natural  to  refer  to  the  tree  structures  which  must  be  al- 
tered at  execution  time  as  dynamic,  and  those  which  remain  struc- 
turally invariant   as  static.    Whether  a  tree  structure  is 
dynamic  or  static  is  independent  of  whether  or  not  the  partitioning 
is  dynamic  or  static.   In  fact,  most  dynamic  partitioning  can  be 
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realized  with  a  static  tree  structure. 

In  order  to  illustrate  some  of  the  ideas  discussed  so  far, 
we  construct  an  example  using  several  levels  of  static  partition 
ing .   Figure  2a  shows  the  partitioning  of  the  array  A  which 
corresponds  to  the  statements  below,  while  Figure  2b  shows  a 
representation  of  the  tree  structure. 

LET  A  BE  A  LOWER  TRIANGULAR  MATRIX  OF  ORDER  (15); 
PARTITION  A  AFTER  ROWS  6,9;  SET  B  =  A<2,1>,  C  =  A<2,2>, 
AND  D  =  A<3,2>;  SET  U  TO  THE  STRICTLY  UPPER  TRIANGULAR 
PART  OF  A<3,1>; 

PARTITION  A<3,1>  AFTER  ROW  1  AND  COLUMN  1;  SET  L  TO 
THE  LOWER  TRIANGULAR  PART  OF  A<3,1><2,2>,  E  TO  THE 
DIAGONAL  PART  OF  A<3,1><2,2>; 

PARTITION  E  AFTER  ROW  2;  SET  F  TO  E<2,2>; 

The  first  partition  statement  divides  the  original  array  A 
into  six  subarrays  which  can  be  referenced.   The  subarrays  are 
subscripted  as  though  A  were  square  and  not  lower  triangular. 
Therefore,  A<1,2>,  A<1,3>,  and  A<2,3>  are  undefined  and  cannot 
be  referenced.   This  subscripting  convention  applies  to  all  of 
the  arrays  in  Table  I.   Also  observe  that  the  subarrays  on  the 
diagonal  inherit  the  lower  triangular  attribute  from  the  parent 
array  A. 

The  subarray  A<3,1>  is  not  specified  in  a  SET  statement; 
however,  its  subscripts  are  constants  so  it  is  equivalent  to  an 
explicit  specification.   For  this  reason  the  corresponding  tree 
structure  is  static. 
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FIG.  2.   Multilevel  partitioning  of  a  lower  triangular 
array:   (a)  array  and  its  partitioned  subarrays ;  (b)  cor- 
responding tree  structure.   The  nodes  for  A<3,1>  and 
A<3,1><2,2>  are  not  labeled  since  they  are  not  explicitly 
renamed. 
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The  notation  A<3,1><2,2>  in  the  second  SET  statement 
specifies  the  subarray  whose  position  is  defined  by  the  sub- 
scripts <2,2>  in  the  parent  array  A<  3 ,  1>  .   This  is  illustrated 
in  Figure  2a. 

Array  Control  Block 

Each  node  of  a  tree  structure  is  referred  to  as  an  array 
control  block  (ACB) .   The  purpose  of  this  control  block  is  to  hold 
control  information  for  accessing  the  elements  of  an  array.   In 
general,  we  are  interested  in  accessing  all  elements  of  a  prescribe 
array  or  subarray;  we  consider  this  to  be  an  intrinsic  property  of 
array  languages. 

Figure  3  shows  an  ACB  node  with  its  fields,  while  Table  II 
lists  some  of  their  characteristics.  In  general,  the  size  of  the 
ACB  node  depends  upon  the  dimension  of  the  array.   Since  we  have 
defined  partitioning  only  for  one  and  two  dimensional  arrays  we 
shall  limit  our  discussion  to  these  cases.   However,  it  should  be 
noted  that  higher  dimensional  arrays  can  be  represented  by  a  simila 
ACB  structure,  and  that  the  corresponding  tree  structures  consist 
of  only  the  root  nodes.   Furthermore,  vectors  (column  or  row)  are 
special  cases  of  rectangular  arrays,  and  need  not  be  distinguished 
from  more  general  rectangular  arrays,  except  by  the  value  of  the 
dimension  field.   The  same  argument  may  be  applied  to  scalars 
which  result  from  partitioning,  since  they  are  also  special  cases 
of  rectangular  arrays.  A  scalar  which  is  not  simultaneously  an 
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FIG.  3.   Array  control  block  (ACB)  node 


TABLE  II   Array  Control  Block  (ACB)  Fields 
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SYMBOL 

FIELD 

ATTRIBUTE 

n 

Name 

Character  String 

A 

Dimension 

Integer 

T 

Type 

Integer 

6D 

Diagonal  Increment 

Integer 

6R 

Row  Increment 

Integer 

CO 

Origin 

Pointer 

Xi 

Lower  Bounds 

Integer 

yi 

Upper  Bounds 

Integer 

P 

Partition 

Pointer 

15 
element  of  an  array  is  not  described  by  an  ACB,  but  is  treated  in 
a  normal  manner.   Therefore,  we  focus  our  attention  upon  two  di- 
mensional arrays  only. 

The  root  node  is  formed  from  information  that  is  contained 
in  a  LET  statement;  in  fact,  we  may  set  up  a  one  to  one  corres- 
pondence between  the  explicitly  declared  arrays  in  the  LET  state- 
ment and  the  root  nodes  of  the  tree  structures.   To  simplify  the 
discussion,  we  shall  assume  that  the  root  nodes  are  formed  upon 
block  entry  and  are  destroyed,  along  with  their  subtrees,  only 
upon  block  exit.   This  condition  will  be  relaxed  later. 

Most  of  the  information  for  filling  a  root  node  is  easily 
obtained.   For  instance,  the  name  (n) ,  dimension  (A),  type  (x)  , 
lower  bounds  (X.),  and  upper  bounds  (y . )  are  explicitly  specified 
in  a  LET  statement,  or  are  implicitly  specified  as  in  the  case  of 
the  type  field  where  a  default  of  rectangular  is  assumed.   The 
partition  pointer  p,  which  is  not  used  until  partitioning  is 
attempted,  is  set  to  the  null  pointer  A. 

In  order  to  define  the  origin,  w,  we  adopt  the  convention  that 
X-=l  for  each  dimension  and  for  all  arrays  and  subarrays.   Then  we 
define  w  as  the  location  of  the  first  accessible  element  of  the 
array  described  by  the  ACB  node.   That  is,  beginning  a  labeling 
of  any  two  dimensional  array  or  any  subarray  with  row  1  and  column  1, 
we  let  w  point  to  the  element  in  position  (1,1)  unless  its  type  is 
strictly  lower  triangular  or  strictly  upper  triangular,  in  which 
case  we  let  it  point  to  the  element  in  position  (2,1)  or  (1,2) 
respectively . 
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TABLE  III   Initial  Values  for  6D  and  6R  in  Root  Nodes 


ARRAY  TYW 

iSTMBOL 

6D 

6R 

Strictly  Lower 
Triangular 

(SLT) 

1 

0 

Lower 

Triangular 

(LT) 

1 

1 

Diagonal 

(D) 

0 

0 

Strictly  Upper 
Triangular 

(SUT) 

-1 

y2-x2-i 

Upper 

Triangular 

CUT) 

-1 

y2-A2 

Tridiagonal 

CTD) 

0 

2 

Rectangular 

(R) 

0 

y2-X2+l 
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Finally,  for  the  root  nodes  only  we  define  the  ACB  fields 
6R  and  6D  from  the  entries  in  Table  III.   We  refer  to  co ,  6R, 
and  6D  as  control  fields,  and  we  consider  their  derivation  in 
the  next  section. 

The  Control  Fields  6R,  6D,  and  co 

The  object  of  this  section  is  to  show  that  the  elements 
of  any  array  or  subarray  may  be  accessed  efficiently  by  using 
the  information  in  the  ACB  nodes .   We  shall  distinguish  two 
cases:   access  of  contiguous  elements  and  access  of  an  ar- 
bitrary element.   For  the  contiguous  case,  the  row  increment, 
6R,  and  the  diagonal  increment,  6D,  are  used  as  initial  values 

We  define  6D  in  Table  III  as  the  number  of  elements  which 
can  be  accessed  in  row  i+1  minus  the  number  of  elements  which 
can  be  accessed  in  row  i.   From  Table  I  it  follows  that  6D 
is  constant  for  all  i,  and  can  be  determined  at  compile  time 
from  the  array  type . 

In  order  to  define  6R,  we  introduce  the  following  nota- 
tion.  Let  S(t)  be  the  set  of  subscript  pairs  (i,j)  which  are 
defined  by  the  relations  in  the  last  column  of  Table  I ,  where 
x  is  the  array  type.   Then  (i,j)  e  S(t)  means  that  the  array 
element  in  position  (i,j)  is  defined  and  can  be  accessed. 
If  (i,j)  £  S(t)  for  some  array  (or  subarray),  then  the  loca- 
tion of  this  element  is  denoted  by  L(i,j). 
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Under  the  conventional  row  major  order  storage  scheme  and 
for  root  nodes  only,  we  define  <5R  in  Table  III  by 

6R  =  L(2,j)  -  L(l,j)         (2,j),  (l,j)  e   S(x) 

for  all  x,  except  diagonal  and  strictly  lower  triangular,  which 
are  defined  as  zero.   It  is  easy  to  see  from  Table  III  that  6R 
can  be  constructed  from  the  information  that  can  be  extracted 
from  a  LET  statement. 

As  indicated  previously,  go,  6R,  and  SB   are  constructed  so 
as  to  access  successive  elements  of  an  array.   It  is  obvious 
that  the  next  element  in  row  i  can  be  accessed  efficiently  by 
means  of 

L(i,j+1)  =  L(i,j)  +  1        (i,j),  (i,j+l)  e    S(x)      (1) 

for  all  t;  this  relation  also  holds  for  any  subarray,  since  the 
row  elements  in  any  subarray  must  necessarily  be  continguous. 
Diagonal  arrays,  of  course,  do  not  have  a  next  row  element  so  (1) 
is  never  satisfied. 

For  accessing  successive  elements  in  any  column,  let 

6R(1)  =  6R,     6R(i+l)  =  6R(i)  +  6D  (2) 

for  all  ACB  nodes.   Then,  the  relation  for  the  location  of  the 
next  element  in  column  j  is  defined  by 

L(i+l,j)  =  L(i,j)  +  6R(i)      (i+l,j),  (i,j)  e  S(x)       (3) 
for  all  x . 
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Formulas  (2)  and  (3)  are  valid  not  only  for  arrays  described 
by  root  nodes,  but  for  all  of  their  subarrays  as  well.   The  veri- 
fication for  arrays  described  by  root  nodes  follows  from  the 
values  of  6R  and  6D  given  in  Table  III,  while  the  verification  for 
subarrays  requires  a  proof  which  we  now  present. 

Let  A  be  a  subarray  of  A   such  that  the  node   is  the  parent 
s  '     p  P        F 

of  node  .   By  definition,  we  say  that  A   is  a  subarray  of  A   if 

and  only  if  every  accessible  element  of  A  ,  considered  as  a  set, 

is  also  an  accessible  element  of  A  ,  in  other  words  A  <=-   A  .   This 

p»  s  -   p 

excludes,  for  example,  A  being  the  tridiagonal  part  of  an  upper 

triangular  A  because  the  elements  on  the  lower  diagonal  of  A  are 
p  s 

not  contained  in  A  . 

P 

We  use  induction  on  the  tree  level,  where  node   is  at  level 

P 

k  and  node   is  at  level  k+1 .   The  case  k=l,  the  root  node  level, 

follows  from  Table  III.   Thus,  assume  that  (2)  and  (3)  are  valid 

for  node  .   Then,  it  follows  that 
P 

6RS  =  5Rp  +  (i-l)6Dp  (4) 

and 


6DS  =  6Dp  (5) 

where  i  is  the  i-th  row  of  A   and  simultaneously  the  first  row 
of  A  .   By  definition  of  a  subarray,  i  always  exists  and  is  equal 
to  or  greater  than  unity.   From  (2)  and  (4)  it  follows  that 

6R  fi)  =  6R   +  (i-l)6D   =  6R  (6) 

ppps  K    J 
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Now  choose  an  arbitrary  element  of  A  and  denote  its  sub- 
scripts by  (m,j).   We  will  show  that  (3)  is  valid  for  A  provided 
that  (m,j)  and  (m+l,j)  e  S  (t) ,  where  S  (x)  corresponds  to  the  set 
of  subarray  elements  of  A  which  can  be  accessed.   Since  the  ar- 
bitrary element  of  A   is  contained  in  A  ,  its  subscripts  with 

s  p 

respect  to  A  are  (i+m,j')  e  S  (t)  for  some  j ' .   Therefore,  since 
(3)  and  (2)  are  valid  for  node   it  follows  from  (6) ,  (5) ,  and  (2)  t 

Ls(m+l,j)  =  Lp(i+m+l,j')  =  Lp(i+m,j')  +  6Rp(i+m) 

=  Ls(m,j)  +  5Rp(i)  +  (m-l)6Dp 

=  Ls(m,j)  +  6Rs  +  (m-l)6Ds 

=  Ls(m,j)  +  6Rs(m). 

Since  the  element  in  A  was  arbitrary,  (3)  is  valid 

for  accessing  any  consecutive  column  elements  of  A  ,  and  this 

completes  the  proof.   This  result  is  summarized  in  Theorem  1. 

We  note  that  (4)  and  (5)  define  the  rules  for  constructing 
the  quantities  6R  and  6D  for  any  subtree  node.   The  only  remaining 
unknown  is  the  value  of  i,  but  this  is  obtained  from  the  parti- 
tioning information  that  is  discussed  in  the  next  section. 

So  far  we  have  shown  that  successive  elements  in  any  row  or 
any  column  can  be  accessed  by  using  the  information  in  an  ACB 
node,  along  with  (1)  or  (3).   Many  other  alternatives  also  exist. 
For  example,  the  successive  elements  along  a  diagonal  are 
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accessed  by  using 

L(i+l,j+l)  =  L(i,j)  +  6R(i)  +  1  (7) 

provided  (i+l,j+l)  and  (i,j)  e  S(x). 

Next  we  consider  how  the  origin  of  a  subarray  is  determined 
from  information  contained  in  the  parent  node.   As  before  we 
shall  qualify  the  arrays  and  fields  by  s  and  p. 

Depending  upon  the  array  type,  we  have  by  definition 

wp  =  Lp(i,j)  (8) 

where  (i,j)  e  {  (1 ,1)  ,  (1 ,  2)  ,  (2  , 1) } .   Since  A   is  a  subarray  of  A  , 
there  exist  integers  6i>_0   and  6j>0  such  that 

ws    =   Lp(i  +  <SiJ+<SJ).  (9) 

Therefore,    it   follows    from    (1)    and    (3)    that 

Lp(i+6i,j+6j)    =    Lp(i+6i,j)    +    5j 

i+6i-l 
■    L    (i,j)    +        I        <5R    (k)    +    6j. 
^  k  =  i  p 

From  the  definition  of  6R(i)  it  follows  that  for  any  node 

i+6i-l 

I        6R(k)  =  6R(i)6i  +  6D(6i  (6i -1) )/2  (10) 

k=i 

where  i=l  or  2,  depending  upon  the  type  t.   Therefore, 

^s  =  wp  +  ^pCi)61  +  6D  (6i(6i-l))/2  +  6j  (11) 

for  all  t. 
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The  unknown  quanities  6i  and  6j  are  determined  from  the 
partitioning  information  which  is  discussed  in  the  next  section. 
Since  by  definition,  i  =  1  or  2 ,  we  may  replace  6R(i)  with  either 
6R(1)  =  6R  or  6R(2)  =  6R  +  6D.   Therefore  (8)  shows  that  the 
origin  of  any  subarray  co   can  be  determined  from  the  information 
in  the  parent  node  and  from  appropriate  partitioning  information. 
We  summarize  these  results  in  the  following  theorems: 

Theorem  1.   (Contiguous  Access)    If  the  fields  6R  and  6D 
of  a  root  node  of  a  tree  structure  are  filled  according  to 
Table  III,  and  if  the  fields  of  each  subtree  node  are  derived 
from  the  parent  node  by 

6RS  =  6Rp  +  (i-l)6D  (12) 

6Ds  =  6Dp  (13) 

where  i  is  the  i-th  row  of  the  array  described  by  node   and  also 
the  1-st  row  of  the  array  described  by  node  ,  then  we  may  locate 
the  successor  of  any  array  element  described  by  a  node  in  any  row, 
column,  or  diagonal  by  applying  respectively 

L(i,j+1)  =  L(i,j)  +  1,  (i,j+l),(i,j)  e  S(x)   (14) 

L(i+l,j)  =  L(i,j)  +  6R(i),  (i+l,j),(i,j)  e  S(x)   (15) 

L(i+l,j+l)  =  L(i,j)  +  6R(i)  +  1      (i+l,j+l),(i,j)  e  S(x) 

where  t  is  any  array  type  and   R(i)  is  defined  by  (2) . 
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Theorem  2.   (Arbitrary  Access)    Let  the  conditions  of 
Theorem  1  be  satisfied,  and  let  L(i,j)  be  the  location  of  an 
arbitrary  element  in  an  array  described  by  a  node.   Then,  for 
any  5i^0  and  6j>_0   such  that  (i  +  6i,j+6j)  e  S(t), 

L(i+6i,j+6j)  =  L(i,j)  +  6R(i) 6i+6D(6i (6i- 1) ) /2  +  6j     (17) 

for  all  array  types  t,  where  <5R(i)  is  defined  by  (2). 
It  follows  from  (17)  that 

L(i+6i,j+6j)  =  co   +  6R(i)6i  +  6D(6i  (6i-l)  )/2  +  6j       (18) 

where  co   =  L(i,j)  and  (i,j)  e  {(1,1),  (1,2),  (2,1)}  depending 
upon  the  type  x.   Therefore,  an  arbitrary  element  of  any  array 
can  be  located  by  using  the  control  fields  co,  6R,  6D  along  with 
the  displacements  6i  and  6j  which  can  be  determined  from  the  sub- 
scripts of  a  referenced  element,  say  A(k,m).  (The  round  parentheses 
refer  to  an  element  of  an  array,  whereas  A<k,m>  refer  to  a  sub- 
array  of  A.)  It  also  follows  that  (11)  is  a  special  case  of  (17) . 

Partition  Control  Blocks 

If  an  array  is  described  by  an  ACB  and  is  not  partitioned, 
then  the  partition  field  p  is  set  to  A,  the  null  pointer,  and  the 
ACB  node  becomes  a  terminal  node  in  the  tree.   On  the  other  hand, 
if  it  is  partitioned,  then  p  is  set  to  point  to  a  partition  con- 
trol block  (PCB)  which  contains  the  necessary  information  for 
deriving  the  nodes  for  the  next  level  in  the  tree  structure. 
Structurally,  the  PCB  may  be  considered  as  a  field  in  the  ACB 
with  a  variable  number  of  subfields;  however,  we  prefer  to  treat 
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FIG.  4.   Partition  control  block  (PCB) 
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it  as  a  separate  control  block. 

The  general  structure  of  the  PCB  is  shown  in  Figure  4.   The 
field  v  contains  the  number  of  partition  lines  defined  by  the 
partition  statement,  namely  M+N,  while  the  fields  r  ,  c   and 
rN+l  »  CM+1  contai-n  tne  lower  bounds  and  upper  bounds,  respectively, 
of  the  array  being  partitioned. 

The  partitioning  lines,  as  indicated  earlier,  are  defined 
by  expressions  in  the  partition  statement.   Corresponding  to  each 
row  expression  R-  and  each  column  expression  C- ,  there  is  a  row 
partition  field  r.  and  a  column  partition  field  c.  in  the  PCB. 
Implicit  in  the  number  of  partition  expressions  is  also  the  maxi- 
mum number  of  subarrays  resulting  from  the  partitioning  of  the 
array.   Therefore,  at  compile  time  one  has  the  information  that 
is  necessary  for  constructing  all  of  the  subtree  nodes  corresponding 
to  a  given  partition  statement.   However,  under  the  assumption 
that  only  a  few  of  these  subarrays  will  be  referenced  by  the 
user,  we  defer  the  building  of  the  subtree  nodes  until  directed 
to  do  so  by  an  explicit  specification,  such  as  "SET  C=A<2,1>." 

The  subarray  pointer  list,  a,  in  the  PCB  is  used  to  hold 
the  pointers  which  point  to  the  subtree  nodes  at  the  next  level 
of  the  tree  structure.   We  can  now  visualize  the  basic  steps  in 
deriving  the  information  for  a  subtree  node.   For  example,  suppose 
that  the  subarray  A<2,1>  is  specified  in  a  SET  statement  and  a 
PCB  for  A  exists.   Then  the  subarray  subscripts  <2,1>  are  used 
to  uniquely  reference  a  pointer  in  a,  which  points  to  the  ACB 
node  for  A<2,1>.   Next,  the  subscripts  <2,1>  are  used  to  uniquely 
reference  the  partition  lines  r,  ,  r2,  and  cn,  c-,  in  the  PCB, 
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which  define  the  boundaries  of  A<2,1>.  With  this  information 
along  with  that  from  the  ACB  node  for  A  we  can  fill  the  ACB  node 
for  A<2,1>.   A  more  precise  description  will  be  presented  in  the 
section  on  basic  algorithms. 

The  pointer  list  tt  in  the  PCB  is  used  to  point  to  the  ACB 
subtree  nodes  which  describe  some  "part  of"  an  array,  such  as 
the  upper  triangular  part  of  A.   For  simplicity,  we  take  it  to  be 
an  array  of  six  pointers;  one  for  each  of  the  nonrectangular  types 
in  Table  I.   We  use  the  notation  tt(t)  to  refer  to  the  pointer 
field  which  corresponds  to  an  array  of  type  t.   If  a  particular 
part  of  an  array  is  not  desired,  then  the  corresponding  pointer 
tt(t)  is  set  to  A.   If  diagonal  partitioning  is  not  desired,  then 
the  pointer  list  tt  need  not  exist.   On  the  other  hand,  if  only 
diagonal  partitioning  is  desired,  then  the  row  partition  fields 
r.,  column  partition  fields  c.,  and  the  pointer  list  a  need  not 
exist . 

The  previous  discussion  illustrates  the  compactness  of  the 
partitioning  information  in  a  PCB,  and  it  also  illustrates  the  gen- 
eral technique  for  adding  subtree  nodes.   Furthermore,  it  should  be 
clear  that  the  PCB  contains  sufficient  partitioning  information  for 
dynamically  adding  a  node  to  the  tree  structure  at  execution  time. 
This  happens,  for  instance,  if  a  subarray  reference  of  the  form 
A<I,J>  is  encountered,  where  I  and  J  are  variables.   In  this  case, 
separate  ACB  nodes  are  created  for  each  distinct  subscript  pair 
<I,J>  as  they  are  needed.   It  is,  of  course,  possible  to  use  a 
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FIG.    6.      Multilevel   partitioning   of  a    lower   triangular   array 
and   relative    addresses    of   the   elements. 
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single  ACB  node  for  all  A<I,J>,  but  we  have  not  chosen  this 
alternative . 

Figure  5  illustrates  the  tree  structure  corresponding  to 
Figure  2  using  simplified  ACB  nodes  and  partition  control  blocks. 
Most  of  the  information  displayed  in  Figure  5  follows  from 
Figure  6  and  from  the  example  which  defined  Figure  2.   A  more 
complete  description,  including  a  derivation  of  the  fields  not 
shown,  is  discussed  in  the  next  section. 

Basic  Algorithms 

The  purpose  of  this  section  is  to  define  when  an  ACB  node 
should  be  added  to  a  tree  structure  and  how  the  information  for 
a  new  node  is  to  be  derived.   We  define  these  processes  in  terms 
of  basic  algorithms. 

In  order  to  interpret  the  examples  presented  in  this  section 
and  in  the  appendix,  we  define  the  following  semantic  rules   which, 
in  general,  are  not  included  in  the  basic  algorithms. 

1.  LET  statements  declare  arrays  and  therefore  are  non- 
executable.  Within  a  block  all  LET  statements  are  processed 
before  any  executable  statement. 

2.  PARTITION  and  SET  statements  are  executable. 

3.  Once  a  PARTITION  statement  is  executed,  any  partitioned 
array  in  that  statement  remains  partitioned  until  control  passes 
out  of  the  block,  or  until  another  partition  statement  in  the 
same  block  redefines  the  partitioning  for  that  array  (an  excep- 
tion is  explained  in  5).   If  a  new  partitioning  is  defined  for 
that  same  array,  then  the  subtree  for  the  previous  partition  is 
deleted  and  a  new  subtree  is  built  using  the  basic  algorithms. 
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If  the  same  partition  statement  is  executed  more  than  once  and  no 
other  partition  statement  redefines  the  partitioning  for  that 
same  array,  then  the  subtree  remains  intact  and  the  fields  of 
the  specified  ACB  or  PCB  are  merely  updated. 

4.  Once  a  SET  statement  is  executed  such  as,  "SET  X=A<2,1>," 
the  subarray  name  X  may  be  used  until  control  passes  out  of  the 
block,  or  until  another  SET  statement  uses  the  same  name  for  a 
different  subarray  (an  exception  is  explained  in  5) .   If  the  name 
X  is  "switched"  to  another  subarray,  then  the  previous  subarray 
A<2,1>  cannot  be  referenced,  but  its  ACB  node  is  not  destroyed. 
This  allows  the  same  name  X  to  be  used  for  different  subarrays 
without  the  overhead  associated  with  creating  and  destroying 
nodes  and  subtrees.   If  the  same  SET  statement  is  executed  more 
than  once  and  no  other  SET  statement  uses  the  same  name,  then 
subsequent  executions  will  update  the  subtree  node  for  X.   How- 
ever, if  the  subarray  is  A<I,J>,  then  a  new  node  will  be  created 
whenever  the  values  of  the  subscripts  <I,J>  reference  a  new  sub- 
array  A<I,J>.   In  this  case  the  name  X  always  refers  to  the  last 
subarray  referenced.   Thus,  a  single  SET  statement  which  is 
executed  repetitively  can  generate  many  subtree  nodes,  one  for 
each  distinct  pair  of  values  of  <I,J>. 

5.  A  SET  statement  of  the  form  SET  A=B ,  where  B  is  defined 
by  some  previous  statement,  is  treated  in  a  different  manner  than 
described  above.   In  particular,  an  ACB  node  is  created  for  A, 
which  is  an  exact  copy  of  the  ACB  node  for  B.   Then,  the  ACB  node 
for  A  becomes  a  root  node  of  a  new  tree  structure,  its  "descendent 
tree,"  which  is  independent  of  the  tree  structure  which  contains 
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the  ACB  node  for  B.   If  another  SET  statement  redefines  A,  then 
the  existing  ACB  node  for  A  and  its  subtree  nodes  are  deleted  and 
a  new  node  is  created  using  the  basic  algorithms.   On  the  other 
hand,  if  the  same  SET  statement  is  executed  more  than  once,  then 
the  ACB  for  A  is  updated,  provided  A  has  not  been  used  in  another 
SET  statement  in  which  case  an  ACB  node  for  A  must  be  recreated 
using  the  basic  algorithms. 

The  notation  which  is  used  in  the  following  algorithms  is 
similar  to  that  found  in  Knuth  [6],  however,  at  times  it  is 
necessary  to  introduce  additional  notation.   This  is  explained 
either  within  the  algorithm  or  in  the  discussion  which  follows 
it. 


32 

ALGORITHM  R   (Root  Node)    The  root  node  of  a  tree  structure 
is  constructed  in  this  algorithm.   It  is  a  straightforward  initial- 
ization of  the  root  node  fields  from  information  contained  in  a 
LET  statement  as  discussed  in  previous  sections. 

We  denote  by  A  the  name  of  the  array  and  by  A"  the  explicit 
pointer  variable  that  is  used  to  reference  the  ACB  which  des- 
cribes A. 

Rl.   [Initialize]   Allocate  ACB,  set  R.   The  size  of  the  ACB  is 
actually  a  function  of  the  dimension  of  the  array,  which  is 
available  at  compile  time.   Once  allocated  the  pointer  R  is 
set  to  point  to  this  ACB  node,  which  is  now  designated  by 
node (R) .   Similarly,  the  fields  of  node(R)  are  designated 
by  n(R),  A(R),  t(R),  ...,  p(R),  where  the  symbols  are  defined 
in  Table  II.   Allocate  A,  set  o>(R)  .   This  allocates  memory 
for  the  array  elements ;  the  array  bounds  and  the  array  type 
determine  the  amount  of  memory  required.   w(R)  is  the  loca- 
tion of  the  first  stored  element  of  the  array  (for  the  two 
dimensional  arrays  of  Table  I  this  corresponds  to  either  the 
(1,1),  (1,2),  or  (2,1)  element). 
R2 .   [Fill  ACB]   n  (R)  «-  name  of  array,  x  (R)  «-  type  of  array  (de- 
fault is  rectangular),  A(R)  «-  dimension  of  array,  X.(R)  «-  1 
lower  bound  for  i-th  dimension  (defined  to  be  unity  for  each 
dimension),  li-(R)  «-  upper  bound  for  i-th  dimension,  p(R)  ■*-   A 
(null  pointer).   Set  <5R(R)  and  6D(R)  as  defined  in  Table  III. 
R3 .   [Set  ACB  pointer]   A«-R,  terminate  algorithm. 
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If  other  arrays  are  declared  in  the  same  LET  statement  or 
if  there  are  other  LET  statements  in  the  same  block,  then  this 
algorithm  is  repeated  until  all  root  nodes  are  filled.   Once 
this  is  completed  the  remaining  algorithms  can  be  used  to  generate 
subtree  nodes  for  the  tree  structures  as  needed. 

As  an  example  consider  the  LET  statement  which  generates 
the  root  node  in  Figures  2  and  5,  namely 

LET  A  BE  A  LOWER  TRIANGULAR  MATRIX  OR  ORDER  (15). 
Using  the  relative  addressing  indicated  in  Figure  6  and  applying 
algorithm  R,  we  obtain  the  first  line  of  Table  V.   The  other 
entries  in  Table  V  are  the  result  of  applying  subsequent  algorithms 
and  are  discussed  at  the  appropriate  place. 
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ALGORITHM  P   (Partition  Array)    Let  P  be  a  pointer  to  an 
ACB  node  in  a  tree  structure,  and  let  a  PARTITION  statement  be 
defined  for  the  array  described  by  node(P).   Denote  by  RT  (CKI<N) 
and  C,  (0<J<_M)  the  row  and  column  expressions  in  the  PARTITION 
statement,  and  assume  that  either  N^O  or  M^O  so  that  we  have  at 
least  one  partition  line.   Then,  this  algorithm  fills  the  PCB 
with  the  necessary  partitioning  information.   If  the  PCB  doesn't 
exist,  it  is  created  in  step  PI.   Q  is  a  pointer  which  is  used  to 
reference  the  PCB. 

PI.   [Locate  PCB]   If  t(P)^R  then  N«-max(N,M),  M^N.   Set  Q«-p(P). 
If  Qj^A  then  go  to  step  P2;  otherwise  allocate  PCB,  set  Q, 
p(P)«-Q,  tt(Q)<-A,  and  a(Q)«-A.   (The  number  of  pointers  in 
a(Q)  depends  on  the  type  t (P)  and  the  number  of  partition 
lines,  namely  M+N) .   Set  v(Q)=M+N. 
P2.   [Repartition?]   If  v(Q)=M+N  then  go  to  step  P3;  otherwise 
delete  subtrees  of  node(P)  which  are  pointed  to  by  a(Q). 
Save  tt(Q),  reallocate  other  parts  of  PCB,  v(Q)^M+N,  a(Q)«-A. 
(The  size  of  the  PCB  is  adjusted  when  repartitioning  of  an 
array  occurs  but  tt(Q)  is  left  unaltered.) 
P3.   [Fill  PCB]   If  x(P)  is  not  rectangular  in  step  PI,  either 
N  or  M  is  zero,  since  either  RT  or  Cj  is  specified  but  not 
both.   In  such  a  case,  interpret  RK-=CK  for  l<K<max  (M,N)  in 
the  following  assignments.   Set  r  (Q,  O)"*-^- (P) -1 ,  r(Q,I)«-Rj 
1<I<N,  r(Q,N+l)^y1(P)  .   Similarly,  c  (Q,0)«-A2(P) -1 ,  c(Q,J)^Cj 
1<J<M,  c(Q,M+l)^-u2(P)  . 
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P4.   [Order  partitions]   Sort  r(Q,I)  1^0,  I^N+1,  and  c(Q,J)  J^O, 

J^M+1  into  ascending  order.   At  execution  time,  the  partition 
lines  are  ordered  according  to  the  values  of  the  expressions 
RT  and  C,  and  not  according  to  the  subscripts  I  and  J. 
Terminate  algorithm. 

If  other  arrays  are  to  be  partitioned  with  the  same  partition 
statement,  then  this  algorithm  is  repeated  the  necessary  number 
of  times.   In  each  case  we  assume  the  existence  of  an  algorithm 
which,  if  necessary,  traverses  the  tree  at  execution  time  to  find 
the  correct  ACB  node  and  set  the  pointer  P.   For  instance,  if  A 
is  the  name  of  the  array  to  be  partitioned,  then  a  tree  traversal 
is  not  necessary  because  the  explicit  pointer  A"  exists  and  P-«-A 
locates  the  correct  ACB  node.   On  the  other  hand,  if  A<2,1>  is 
the  array  to  be  partitioned,  then  no  explicit  pointer  exists  for 
the  ACB  node  for  A<2,1>.   Therefore,  we  traverse  a  subtree  by 
starting  with  A,  then  locating  the  PCB  through  p (A) ,  and  finally 
locating  the  ACB  for  A<2,1>  by  using  cr(Q,2,l)  which  is  the  unique 
pointer  field  in  a(Q)  that  is  reserved  for  the  ACB  for  A<2,1>.   The 
pointer  Q  is,  of  course,  equal  to  p (A) .   Thus,  in  this  case  P 
would  assume  the  value  of  a(Q,2,l).   This  general  procedure  can 
be  extended  to  more  than  two  levels.   Thus,  whenever  an  algorithm 
assumes  that  P  points  to  an  ACB  node,  we  will  assume  that  its 
value  was  obtained  by  a  simple  traversal  algorithm  of  this  type. 

As  indicated,  we  use  a(Q,I,J)  to  refer  to  the  pointer  in 
a(Q)  which  points  to  the  ACB  node  for  A<I,J>,  where  A  is  the  array 
that  is  partitioned.   For  our  purposes  we  may  think  of  cr(Q),  as  a 
two  dimensional  array  of  pointers,  one  for  each  possible  subarray. 
However,  in  practice  the  size  of  a(Q)  is  reduced  to  the  maximum 


arrays  defined  by  N  row  partitions  R,  ,  ... 
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number  of  subarrays  that  can  legally  be  referenced.   For  instance 

if  A  is  lower  triangular  then  there  are  only  (N+l)(N+2)/2  sub- 

Ri 

Consider  now  the  situation  where  another  partition  statement 
repartitions  the  same  array.   The  old  partition  must  be  deleted, 
along  with  all  of  the  subtrees,  and  the  size  of  the  PCB  must  be 
adjusted  to  account  for  a  different  number  of  subarrays.   This 
is  accomplished  in  step  P2. 

If  x  is  lower  triangular,  we  denote  by  tt(Q,t)  the  pointer 
in  tt(Q)  which  points  to  the  ACB  node  which  describes  the  lower 
triangular  part  of  A.   Thus  in  step  PI,  tt(Q)-<-A  simply  initializes 
all  the  pointers  in  tt  (Q)  .   The  pointer  fields  in  tt(Q),  therefore 
are  available  later  for  diagonal  partitioning  (Algorithm  D)  if 
needed. 

Algorithm  P  may  create  a  PCB;  however,  if  diagonal  partition- 
ing is  encountered  first,  then  Algorithm  D  creates  the  PCB.   In 
this  case  o (Q)    is  empty,  (v(Q)=0),  thus  if  a  partitioning  state- 
ment is  encountered  and  a  PCB  exists,  it  becomes  necessary  to 
expand  the  PCB  in  step  P2  without  altering  the  pointers  in  tt(Q). 

Rectangular  arrays  may  be  partitioned  after  rows  and  columns, 
or  after  rows  only  or  after  columns  only.   In  the  case  of  vectors, 
only  the  latter  are  meaningful.   In  such  cases,  we  define 
A<I,0>=A<I>  if  M=0  and  A<0,J>=A<J>  if  N=0.   It  is  also  obvious 
that  in  such  cases  a (Q , I , 0) =a (Q, I)  becomes  a  simple  array  of 
pointers  of  length  N+l  or  M+l.   Partitions  of  this  type  occur  in 
the  appendix,  and  the  algorithms  defined  here  apply  using  the 
above  notational  changes. 
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To  illustrate  Algorithm  P,  consider  the  statement 
PARTITION  A  AFTER  ROWS  6,9 
which  corresponds  to  the  PCB  at  level  1  in  Figure  5.   The  PCB  is 
created  in  step  PI  and  filled  in  step  P3.   Since  x(P)  is  lower 
triangular  the  quantities  C,  and  C2  in  step  P3  are  determined 
from  R,  and  R? .   That  is,  by  definition,  the  partitioning  lines 
R,=6  and  R7  =  9  are  reflected  off  of  the  main  diagonal  giving  C-,=6 
and  C2=9.   It  is  also  important  to  note  that  the  pointers  in  a 
and  tt  are  set  to  A  in  step  PI,  and  it  is  the  responsibility  of 
subsequent  algorithms,  namely  S  and  D,  to  insert  pointers  into 
the  appropriate  places  in  a   and  tt  to  create  the  branches  shown 
in  Figure  5. 
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ALGORITHM  S   (Subtree  Node)   Let  P  be  a  pointer  to  an  ACB 
node,  and  let  a  SET  statement  be  defined  which  sets  an  array, 
say  C,  to  the  subarray  described  by  node(P).   For  example, 
"SET  C  =  A<2,1>."   Then,  this  algorithm  defines  the  rules  for 
constructing  a  subtree  node  for  C.   We  assume  that  P  is  deter- 
mined as  in  Algorithm  P.   Throughout  this  algorithm,  P  points 
to  the  parent  node,  Q  to  the  PCB  for  node(P),  and  S  to  the  ACB 
subtree  node.   The  pointer  P  is  obtained  by  a  simple  traversal  of 
the  tree  if  necessary,  but  for  the  above  example  it  would  be  P«-A. 

51.  [Initialize]   Set  Q  *■   p(P).   We  assume  that  a  PCB  has  been 
filled  by  algorithm  P;  otherwise,  the  partitioning  is 
undefined. 

52.  [ACB  exist?]   S  +■   a(Q,I,J)  where  I  and  J  are  extracted  from 
the  subarray  subscript  pair  <I,J>  in  a  reference  of  the  form 
A<I,J>.   If  S=A  then  allocate  ACB,  set  S,  a(Q,I,J)  «-  S, 
p(S)  «-  A. 

53.  [Fill  ACB]   Set  6D(S)  «-  6D(P)  .   If  I=J  then  t(S)  «-  t(P), 
otherwise  t(S)  +-  R.   Set  Ai(S)  «-  1  for  i=l,2.   y-^S)  «-  r(Q,I) 
-r(Q,I-l),  y2(S)  +  c(Q,J)-c(Q,J-l),  6R(S)  <-  6R(P) 

+6D(P)  r(Q,I).   If  t(P)=SLT  then  6i  «-  r(Q,I-l)-l,  6r  «-  6R(P) 
+6D(P);  otherwise  6i  <-  r(Q,I-l),  6r  <-  6R(P)  .   If  t(P)=SUT 
then  6j  +■  c(Q,J-l)-l;  otherwise  6j  +   c(Q,J-l).   Set  co(S)  ♦■  w(I 
+  6r  6i  +  6D(P) (6i(6i-l))/2  +  6 j .   If  the  attribute  "VECTOR" 
follows  the  subarray  reference  then  set  A(S)«-1  and  go  to 
step  S4 .   If  the  attribute  "SCALAR"  follows  the  reference 
then  set  A(S)^-0,  and  go  to  step  S4  ;  otherwise  set  A(S)«-A(P). 
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If  t(S)=SLT  then  u>(S)  *■  w(S)  +  6R(S).   If  t(S)=SUT  then 

u)(S)  *  a)(S)  +  1. 
S4.   [Name  node]   Set  £   «■  S,  n  (S)  +■  'C,  where  C  is  the  name  of 

the  subarray  which  appears  in  the  SET  statement. 

The  first  time  a  SET  statement  is  executed  an  ACB  is  created 
for  the  subtree  node.   Subsequent  executions  of  the  same  SET 
statement  may  or  may  not  create  another  ACB,  it  depends  upon 
whether  I  and  J  have  values  which  are  distinct  from  the  previous 
execution.   However,  in  many  applications  it  is  sufficient  to 
use  constants  in  place  of  I  and  J,  then  the  ACB  node  is  merely 
updated  on  subsequent  executions.   Examples  which  illustrate 
this  point  are  presented  in  the  Appendix. 

We  note  that  the  pointer  C  is  always  set  in  step  S4  to 
point  to  the  current  ACB  node  for  the  subarray  A<I,J>.   Within 
a  loop,  therefore,  C  may  be  used  to  refer  to  different  subarrays 
of  the  form  A<I,J>.   The  generality  of  this  feature  includes  even 
the  possibility  of  varying  the  type  of  the  subarray.   For  example, 
if  A  is  lower  triangular,  then  A<I,J>is  lower  triangular  if  I=J 
and  rectangular  otherwise,  so  that  the  type  is  clearly  a  function 
of  the  subscript  values. 

The  values  of  <5i  and  6j  in  step  S3  are  obtained  from  the 
PCB  fields  r(Q,I-l)  and  c(Q,J-l),  which  represent  the  current 
position  of  the  partitioning  lines.   Also,  recall  that  Algorithm  P 
insures  that  the  partitioning  lines  are  in  ascending  order;  it 
follows,  therefore,  that  y-(S)>_0  in  step  S3.   If  u-(S)=0  for  some 
i,  then  we  define  the  subarray  which  has  both  X.=0  and  u.=0  to  be 
a  "NULL"  array.   This  is  discussed  in  the  next  section.   The  case 
where  r(Q,J)  assumes  a  negative  value  is  not  allowed. 
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We  can  illustrate  algorithm  S  by  considering  the  statement 
SET  B=A<2,1>,  C=A<2,2>,  D=A<3,2> 
which  builds  the  ACB  nodes  for  B,  C,  and  D  in  Figure  5  at  level  2. 
The  values  of  the  fields  are  shown  in  Table  V.   Each  node  requires 
a  separate  application  of  Algorithm  S. 

The  attributes  VECTOR  and  SCALAR  in  step  S3  play  a  special 
role.   In  essence,  they  allow  the  compiler  to  invoke  a  precedence 
relation  when  compiling  array  expressions  so  that  the  number  of 
operations  can  be  reduced.   The  essential  point  is  that  the  attri- 
bute characterizes  the  shape  of  the  dynamic  subarray  as  it  moves 
over  its  parent  array,  and  the  attribute  guarantees  that  the  shape 
will  not  change.   The  examples  in  the  appendix  show   that  array 
algorithms  can  be  constructed  which  have  dynamic  subarrays  of 
constant  shape.   Though  array  evaluation  is  not  central  to  par- 
titioning, it  is  vital  to  array  language  design  and  deserves  some 
mention  at  this  point. 
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ALGORITHM  D   (Diagonal  Partitioning)   Let  P  be  a  pointer  to 
an  ACB  node,  and  let  a  SET  statement  be  defined  which  sets  an 
array,  say  B,  to  a  "part  of"  the  array  described  by  node(P).   For 
example,  "SET  B  TO  THE  LOWER  TRIANGULAR  PART  OF  .  .  ."is  typical. 
This  algorithm  prescribes  the  rules  for  deriving  the  information 
for  the  subtree  node  from  the  parent  node  and  the  given  SET  state- 
ment . 
Dl.   [Locate  PCB]   Set  Q  +■  p(P).   If  Q^A  then  go  to  step  D2 , 

otherwise  allocate  PCB,  set  Q,  p(P)  «-  Q,  t(Q,t)  «-  A  for  all  t. 
Set  v(Q)=0  (indicates  that  a(Q),  r(Q,I),  and  c(Q,J)  are  empty.) 
D2.   [ACB  exist?]   Set  S  «-  tt(Q,t)  where  t  is  the  type  prescribed 
in  the  SET  statement.   If  S^A   then  go  to  step  D3 ,  otherwise 
allocate  ACB,  set  S,  tt(Q,  t)  *■  S,  p(S)  +■  A. 
D3.   [Fill  ACB]   Set  A(S)  <-  A(P),  6D(S)  +-  6D(P),  6R(S)  «-  6R(P)  , 
X.(S)  «-  1  i  =  l,2,  t(S)  *■   type  of  subarray  prescribed  in  SET 
statement.   Set  y.(S),  and  co(S)  according  to  Table  IV. 
D4 .   [Name  node]   Set  n(S)  ■*•  *B'  where  B  is  the  name  of  the  sub- 
array  given  in  the  SET  statement.   Set  B*-<-S  where  B  is  the 
pointer  variable  which  corresponds  to  the  name  B.   Terminate 
algorithm. 
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TABLE  IV  Rules  for  Calculating  y-(S)  and 
oo(S)  in  Algorithm  D 


t(P) 
t(S) 

SLT 

LT 

D 

SUT 

UT 

TD 

R 

SLT 

X 

2,4 

X 

X 

X 

X 

1,4 

LT 

X 

X 

X 

X 

X 

X 

1,5 

D 

X 

2,5 

X 

X 

2,5 

2,5 

1,5 

SUT 

X 

X 

X 

X 

2,3 

X 

1,3 

UT 

X 

X 

X 

X 

X 

X 

1,5 

TD 

X 

X 

X 

X 

X 

X 

1,5 

R 

X 

X 

X 

X 

X 

X 

X 

Key 


Rule 


1   yi(S)  -  min(y1(P),y2(P)) 


2   y.(S)  <-   y.(P) 


3   o)(S)  *■  oj(P)  +  1 


4   oj(S)  -»-  oo(P)  +  6R(P) 


5   w(S)  +■  co(P) 


i  =  l,2 
i=l,2 


X   Not  Allowed 
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The  condition  that  P  points  to  the  parent  node  in  Algorithm 
D  implies  the  existence  of  a  traversal  algorithm.  For  instance, 
consider  the  statement 

SET  U  TO  THE  STRICTLY  UPPER  TRIANGULAR  PART  OF  A<3,1> 
which  generates  two  ACB  nodes  in  Figure  5,  namely  the  node  for 
A<3,1>  and  the  node  for  U.   In  order  to  obtain  P,  we  assume  that 
the  traversal  algorithm  starts  with  the  explicit  pointer  A,  finds 
the  PCB  through  Q«-p(A),  and  discovers  that  6(Q,3,1),  the  pointer 
for  A<3,1>  is  A.   Thus,  it  invokes  algorithm  S  to  create  and  fill 
an  ACB  node  for  A<3,1>;   only  then  does  it  have  the  correct  pointer 
for  P.   Algorithm  D  can  now  be  used  to  create  the  ACB  node  for  U, 
and  incidentally  a  PCB  as  well.   Table  V  shows  the  values  of  the 
fields  for  these  two  ACB  nodes. 

Another  feature  of  the  above  example  is  that  the  PCB  which 
is  created  in  step  Dl  consists  of  only  tt  .  When  the  next  state- 
ment, namely 

PARTITION  A<3,1>  AFTER  ROW  1  AND  COLUMN  1 
is  encountered,  Algorithm  P  must  adjust  the  size  of  the  PCB  to 
accomodate  the  row  and  column  partitions.   Once  this  is  done,  the 
PCB  takes  the  form  as  shown  in  Figure  5. 

Algorithm  D  is  also  applied  in  creating  the  ACB  nodes  for 
A<3,1><2,2>,  L,  and  E  in  Figure  5  and  Table  V.   The  statement 
which  invokes  Algorithm  D  is 

SET  L  TO  THE  LOWER  TRIANGULAR  PART  OF  A<3,1><2,2>, 
E  TO  THE  DIAGONAL  PART  OF  A<3,1><2,2>. 
The  interpretation  is  the  same  as  described  previously. 
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Table  IV  not  only  describes  the  rules  for  determining  the 
upper  bounds  u-  and  the  origin  10  as  a  function  of  the  array  type, 
but  also  specifies  the  types  of  subarrays  that  are  allowed. 
According  to  our  previous  definition,  the  subarray  may  be  of  the 
same  type  as  the  parent  array.   Table  IV  does  not  allow  this. 
The  reason  is  that  under  diagonal  partitioning  you  would  be 
merely  duplicating  the  parent  ACB  node.   Instead,  we  create  a 
new  root  node  using  Algorithm  E  so  that  we  may  have  an  independent 
description  of  any  array  or  subarray.   The  new  root  nodes  are 
different  from  other  root  nodes  because  they  are  not  created  by 
Algorithm  R.   The  new  root  nodes  in  a  sense,  give   rise  to  des- 
cendent  trees  because  they  inherit  their  original  description 
from  some  tree  node;  however,  once  they  are  created  they  may 
grow  independently  of  their  parent  tree.   In  this  way  we  obtain 
independent  partitions  which  have  also  been  referred  to  as  over- 
laying partitions. 
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ALGORITHM  E   (Equivalent  node)   Suppose  a  SET  statement  of 
the  form  "SET  A=B"  is  given,  where  B  denotes  any  array  (or  sub- 
array)  which  has  been  defined  previously.   We  assume  that  B  is  not 
qualified,  that  is,  there  are  no  attributes  or  subarray  subscripts 
attached  to  B.   Secondly,  we  assume  that  the  explicit  pointer  A 
is  either  A,  which  implies  that  A  has  not  been  used;  or  that  A 
points  to  some  subtree  node  which  implies  that  it  is  to  be  re- 
defined by  this  algorithm,  or  that  A  points  to  a  marked  node  which 
implies  that  the  node  was  previously  created  by  this  algorithm. 
We  have  not  previously  indicated  the  existence  of  a  marking  bit 
in  the  ACB  node,  but  we  assume  that  it  now  exists.   This  algorithm 
duplicates  the  information  of  node(B)  in  node (A) . 
El.   [ACB  exist?]   Set  E«-A\   If  E=A  go  to  step  E2.   If  node(E) 

is  marked  go  to  step  E3 ;  otherwise  go  to  step  E2.   (Any  node 
which  is  marked  must  have  been  created  by  this  algorithm  and 
must  be  a  root  node  of  a  descendent  tree.) 
E2.   [Create  ACB  ]   Allocate  ACB,  set  E,  p(E)«-A,  mark  node(E). 
E3.   [Copy  Fields]   Set  t(E)«-t(B~),  6R(E)«-6R(S)  ,  6D  (E)«-6D(5)  , 

Xi(D)-Xi(S),  yi(E)^-yi($),  A(E)-A(3). 
E4 .   [Name  node]   Set  n(E)«-'A',  A*«-E  and  terminate  algorithm. 

If  the  SET  statement  which  invokes  this  algorithm  is  exe- 
cuted more  than  once,  then  subsequent  executions  will  merely  up- 
date the  ACB  node  for  A  with  the  current  information  in  node  (B) . 
There  are  no  restrictions  imposed  on  B;  it  may  be  an  arbitrary 
array  as  long  as  it  is  defined. 
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After  Algorithm  E  is  executed  for  some  SET  statement,  the 
array  A  may  be  partitioned,  and  any  such  partitioning  is  inde- 
pendent of  any  partitioning  for  B.   This  idea  may  be  extended  to 
include  any  number  of  independent  partitions  of  B  by  simply 
adding  more  SET  statements  of  the  form  described. 

It  is  also  possible  to  execute  the  SET  statement,  and  then 
change  B  either  through  another  SET  statement  or  partition  state 
ment.   If  this  is  done  then  it  is  clear  that  A  continues  to  des- 
cribe the  original  B  and  will  continue  to  do  so  as  long  as 
Algorithm  E  is  not  invoked  a  second  time.   The  first  example  in 
the  appendix  illustrates  an  application  of  Algorithm  E  by  using 
the  statement  SET  L  =  A. 

Null  Operands 

Recall  that,  in  general,  the  position  of  a  partitioning 
line  is  defined  by  a  value  of  an  expression  in  a  partitioning 
statement.   This  implies  that  references  to  the  subarrays  can- 
not be  determined  until  execution  time;  it  also  implies  that 
some  of  the  partitioning  lines  may  be  positioned  outside  of 
the  array  boundaries,  or  that  several  partitioning  lines  may 
be  superimposed.   Each  of  these  problems  must  be  solved  in 
order  to  uniquely  define  the  value  of  array  expressions. 

The  first  problem  is  solved  by  ordering  the  contents  of 
the  r.  and  c.  fields  in  the  PCB  into  ascending  order  during 
the  execution  of  the  partition  statement.   This  uniquely  de- 
fines the  boundaries  of  each  subarray,  provided  all  of  the 
partitioning  lines  which  define  the  subarrays  are  distinct 
and  lie  on  the  boundary  or  inside  of  the  parent  array. 
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If  an  array  is  partitioned  and  if  a  partitioning  line 
lies  outside  of  the  array  boundaries,  then  all  subarrays 
which  are  bordered  by  that  line  are  defined  to  have  a  null 
value.   Similarly,  if  two  partitioning  lines  are  superim- 
posed, then  all  subarrays  which  would  normally  lie  between 
these  two  lines,  if  they  were  not  superimposed,  are  also 
defined  to  have  a  null  value.   In  an  array  expression,  an 
operand  which  has  a  null  value  is  defined  to  be  a  null 
operand . 

In  the  OL/2  language,  a  null  operand  in  an  array  express- 
ion obeys  the  same  rules  as  the  additive  identity.   That  is, 
the  null  operand  behaves  like  a  zero  matrix,  a  zero  vector, 
or  a  zero  scalar,  except  in  scalar  division,  in  which  case  it 
is  equivalent  to  no -operation.   The  main  differences  between 
a  null  operand  and  a  zero  operand  is  that  the  zero  operand 
requires  actual  computation  while  the  null  operand  is  combined 
with  other  operands  without  any  computation. 

One  of  the  interesting  by-products  of  defining  null 
operands  is  that  it  simplifies  the  writing  of  array  algorithms 
In  general,  the  simplification  occurs  in  the  initial  and  final 
steps  of  an  algorithm,  which  normally  require  special  consid- 
eration, but  with  null  operands  are  automatically  included 
in  the  general  case.   The  examples  in  the  appendix  illustrate 
this  point. 
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Conclusion 

The  purpose  of  this  paper  was  to  introduce  a  general  par- 
titioning strategy  for  array  languages.   By  defining  differ- 
ent types  of  arrays  and  various  modes  of  partitioning  we  were 
able  to  obtain  an  efficient  and  powerful  means  of  accessing 
subarrays.   Since  partitioning  is  a  fundamental  operation 
that  is  required  in  implementing  many  array  algorithms,  it  is 
of  general  interest  to  the  user  and  the  designer  of  array 
languages . 

Of  prime  interest  to  the  language  designer  are  the  basic 
algorithms  and  the  design  of  the  information  structure  which 
describes  dynamic  subarrays  of  various  shapes  and  sizes. 
Whenever  possible,  we  have  tried  to  make  the  basic  algorithms 
independent  of  the  semantics  associated  with  OL/2,  so  that 
they  may  be  applied  to  the  design  of  other  array  languages. 
It  is  evident  that  other  basic  types  of  arrays  could  be  added 
to  those  we  have  considered;  for  instance,  triangular  arrays 
which  are  obtained  by  rotating  either  lower  or  upper  trian- 
gular arrays  by  90  degrees.   However,  one  should  not  that  the 
basic  array  types  which  we  have  described  are  adequate  for 
handling  the  standard  types  of  sparse  matrices  which  arise  in 
applications  of  partial  differential  equations,  namely,  band 
matrices  and  block  tridiagonal  matrices.   The  technique  is  to 
realize  that  these  types  of  sparse  matrices  can  be  constructed 
by  taking  combinations  of  the  basic  types  of  arrays.   This,  of 
course,  requires  one  to  introduce  additional  data  types,  but 
the  underlying  information  structure  remains  the  same. 
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For  the  user,  partitioning  provides  a  vital  and  conven- 
ient means  for  accessing  parts  of  arrays.   The  result  is  that 
one  may  now  apply  an  array  language  to  the  problem  of  construct- 
ing array  algorithms,  not  only  for  iterative  methods,  but  also 
for  direct  methods.   For  properly  written  array  algorithms, 
the  number  of  operations  for  any  of  the  iterative  or  direct 
methods  is  the  same  as  for  an  element  language,  and  the  amount 
of  memory  space  is  always  optimal  for  each  array.   Only  in  the 
evaluation  of  array  expressions  is  additional  memory  space 
required,  and  this  can  always  be  automatically  generated  by 
the  compiler. 
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APPENDIX 

The  examples  given  here  are  written  in  the  OL/2  language, 
the  purpose  being  to  illustrate  the  practical  aspects  of  par- 
titioning.  In  both  examples  the  number  of  operations  and  the 
amount  of  storage  is  nearly  optimal,  the  difference  being  that 
one  temporary  vector  is  needed  in  the  evaluation  of  some  of 
the  array  expressions. 

There  are  several  variations  of  the  syntax,  some  more 
compact  than  others.   Here  we  present  a  natural  and  descriptive 
form  which  requires  relatively  little  explanation,  yet  is  accept- 
able to  the  compiler. 

CH0LESKYJ1ETH0D:   PROCEDURE (A ,X ,Z ,N )  ; 

LET   A   BE  A  LOWER  TRIANGULAR  MATRIX  OF 
ORDER  (N);    LET   X,   Y,   AND   Z   BE  VECTORS  OF 
ORDER  (N) ;    SET  L=A; 

FOR  K=1,2,...,N;    PARTITION   A,   Y,   AND   Z 
AFTER  ROWS  K-l,  K;    SET  R=A<2,1>  ROW  VECTOR, 
D=A<2,2>  SCALAR,  M=A<3,1>,  C=A<3,2>   COLUMN  VECTOR, 
U=Y<2>,  V=Z<2>,  AND  W=Y<1>; 

D  =  SQRT(  D  -  (RT ,R' )  ); 
C  =  (  C  -  MxR'  )/D; 
U  =  (  V  -  R*W  )/D; 

END; 
BACK_SUBSTITUTION : 

FOR  K=N,  N-l,...,l;    PARTITION   L,   X,   AND 
Y   AFTER  ROWS  K-l,   K;    SET  C=L<3,2>,  D=L<2,2> 
SCALAR,  S=X<2>,  U=Y<2>,  AND  T=X<3>; 

S  =  (  U  -  C'xT  )/D; 

END  CHOLESKY  METHOD; 
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FIG.    7.      Cholesky  method:       (a)    decomposition  and 
forward    substitution;     (b)    back    substitution. 
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The  Cholesky  method  for  solving  a  system  Ax=z  is  preferred 
over  other  direct  methods  when  A  is  dense,  symmetric,  and 
positive  definite.   Since  A  is  symmetric  it  can  be  stored  as  a 
lower  triangular  matrix.   The  first  part  of  the  algorithm  de- 
composes the  matrix  A  into  the  product  L*L'  and  stores  L  into 
A  without  using  additional  storage.   The  forward  substitution 
is  carried  out  at  the  same  time  so  that  y  contains  the  solution 
of  Ly=z.   Finally,  the  back  substitution  solves  L'x=y  for  x. 
Since  A  is  positive  definite,  division  by  zero  is  not  possible; 
however,  an  ON  condition,  as  in  PL/1,  may  be  inserted  if  desired 
Finally,  observe  that  A  and  L  refer  to  the  same  array,  but  they 
are  partitioned  independently  of  each  other. 

The  prime  denotes  the  transpose  of  a  vector  or  matrix. 
Vectors  which  are  declared  in  LET  statements  are  column  vectors. 
The  END  statement  acts  as  a  delimiter  for  procedure  blocks  and 
FOR  statements  as  in  PL/1,  and  it  may  also  serve  as  a  RETURN  as 
in  PL/1.   The  arithmetic  assumes  the  default  attributes  of 
double  precision  floating  point  except  in  the  case  of  the  unde- 
clared variables  K  and  N  which  are  integers. 

All  array  operations  follow  the  usual  rules  of  linear 
algebra,  except  when  an  array  becomes  a  null  operand.   This 
usually  occurs  in  the  initial  and  final  steps  of  an  algorithm 
when  the  partitioning  lines  force  some  of  the  subarrays  to  lie 
outside  of  the  original  array.   For  instance,  when  K=l  in  the 
first  example,  the  subarrays  R,  M,  and  W  in  Figure  7  become 
null  operands,  and  the  array  expressions  degenerate  into 
D  =  SQRT(D),  C  =  (C)/D,  and  U  =  (V)/D  respectively.   Formally, 


54 


o  o  o  o  I  o|  o  o 

O  O  O  O  I  O  I  o  o 

O  O  O  O  |X  I  o  o 

o  o  o  o  o  o  o 

o  o  o  o'o'o  o 


o  o  o  '  o  '  o 
o  o  o  o  I  o 
o  C  o'y|o 
o  o  o  I  o  I  o 
o  o  o  I  o  I  o 
o  o  o  I  o  I  o 


o  o  o  o  o 

o  o  o  o  o 

o  Q    o  o  o 

o  o  o  o  o 

o  o  o  o  o 

o 

o 


o  o  I  o  I  o  o 

(o) 


o  o 

o  o 

o  M  o 

o  o  o 

o  o  o 

o  o  o 

o  o  o 


o 

o 

o 

o 

o 

o 

o 

o 

0 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 
o 
o 

o 
o 

0 

o 
R 
o 

o 
o 
o 

o    o 
jpjo 

o  ■  o 

1        . 

o 
o 
o 

o 

z 

o 

o 
o 
o 

o 
o 
o 

o 
o 
o 

o 
o 
o 

o 

0 

o 

0 

1 

o 

0 

o 

o 

o 

o 

o 

o 

o 

o 

o 

1 

0 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

"o 

o 

o 

o 

o 

0 

o 

o 

o 

o 

o 

o 

I'o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

lo 

0 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

|o 

lo 

o 

o 

0 

o 

o 

o 

(b) 


o\o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o\ 

o^ 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

^o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

^ 

o< 

o 

o 

o 

u 

o 

o 

o 

o 

o 

o 

o 

s? 

o 

o 

o 

o 

o 

o 

0 

o 

o 

o 

o 

o\ 

.o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

> 

>e 

o 

o 

o 

o 

o 

o 

o 

L 

o 

o 

o 

o\o 

o 

o 

o 

o 

o 

o 
o 

o 
o 

o 
o 

o 
o 

o 
o 

o 
o 

0 
0 

o 

o 

\ 

o 

o 

o 
o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o\o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o^ 

(C) 


FIG.  8.   Crout  decomposition:   (a)  first  level 
rectangular  partitioning  of  A;  (b)  second  level 
rectangular  partitioning  of  C,  Y,  and  M;  (c)  first 
level  diagonal  partitioning  of  A. 
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null  operands  obey  the  same  rules  as  the  additive  identity. 

Since  there  is  an  additive  identity  for  vectors,  matrices,  and 
scalars  the  rules  are  well  defined.   The  only  exception  is 
division  by  a  null  operand  which  produces  a  null  result  rather 
than  an  undefined  operation.   Clearly,  null  operands  can  be 
handled  without  actually  carrying  out  the  computations. 

The  second  example  is  the  familiar  Crout  algorithm.   This 
algorithm   requires  two  levels  of  partitioning  to  decompose  a 
matrix  A  into  the  product  L'U.   Theoretically,  L  is  lower  tri- 
angular and  U  is  unit  upper  triangular,  but  in  practice  the 
diagonal  elements  of  U  are  not  stored.   Therefore,  U  is  de- 
clared as  a  strictly  upper  triangular  matrix.     As  the  algorithm 
progresses,  the  elements  of  L  and  U  are  formed  and  stored  in  A. 
Figure  8  illustrates  the  process  for  some  intermediate  value  of  K 

LU_DECOMPOSITION:   PROCEDURES, N)  ; 

LET   A   BE  A  MATRIX  OF  ORDER  (N) ;    SET  L  TO 
THE  LOWER  TRIANGULAR  PART  OF   A,   AND   U   TO  THE 
STRICTLY  UPPER  TRIANGULAR  PART  OF   A; 

FOR  K=l,2,...,Nj   PARTITION   A   AFTER  ROW  K-l 
AND  AFTER  COLUMNS  K-l,   K;    SET  C=A<2,1>,  X=A<1,2>, 
Y=A<2,2>,  B=A<1,3>,  AND  M=A<2,3^;    PARTITION  C,  Y,  M 
AFTER  ROW  lj   SET  R=C<1>,  D=Y<1>  SCALAR,  AND  Z=M<1>; 

Y  =  Y  -  CxX; 

Z  =  (  Z  -  RxB  )/D; 

END  LU  DECOMPOSITION; 
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Though  L  and  U  are  specified,  they  are  not  used  in  this 
example  except  to  clarify  the  algorithm.   However,  another 
segment  could  be  added  which  would  use  L  and  U  to  solve  one 
or  more  systems  of  equations  Ax-=z.  in  the  usual  way.   Finally, 
partial  pivoting  may  be  introduced  by  using  a  slightly  different 
partitioning  and  an  interchange  statement;  however,  this  does 
not  add  anything  new  to  the  example. 
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